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Glass-forming liquids have been extensively studied in recent decades, but there is still no theory 
that fully describes these systems, and the diversity of treatments is in itself a barrier to understand- 
ing. Here we introduce a new simple model that (possessing both liquid-crystal and glass transition) 
unifies different approaches, producing most of the phenomena associated with real glasses, with- 
out loss of the simplicity that theorists require. Within the model we calculate energy relaxation, 
non- exponential slowing phenomena, the Kauzmann temperature and other classical signatures. 
Moreover, the model reproduces a subdiffusivc exponent observed in experiments of dense systems. 
The simplicity of the model allows us to identify the microscopic origin of glassification, leaving 
open the possibility for theorists to make further progress. 



* This work was carried out while DC was at the School of Chemistry and Chemical Biology, University College Dublin, Belfield, Dublin 
4, Ireland. 
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I. INTRODUCTION 

Understanding the nature of the glass transition is considered to be one of the great outstanding problems of the 
condensed state of matter [1, 2]. Advances have been real, but progressive, and there is still no theory of the glassy 
state that explains all of its properties. Quite different approaches and levels of description have arisen within the 
field, each with acknowledged strengths. However, the diversity of levels of description has in itself become a barrier 
to progress and understanding, for in many cases there is no apparent connection between the different strands of 
research. 

Here we cannot summarize the numerous efforts, but mention only several trends in the field that are relevant to 
this paper. The glass problem is related to a wider set of phenomena, collectively named dynamical arrest: that 
process in which many particles dramatically slow in a concerted manner [3-15]. One interpretation is as follows: for 
simple repulsive interactions, with increasing density, progressive loss of space around a typical particle leads it to 
become effectively trapped by its neighbors, a phenomenon often termed caging [16, 17]. Occasionally, the particle 
can escape from the cage and make longer movements before being trapped in another cage. Kinctically constrained 
models, such as the one introduced by Kob and Andersen [3], represent the intra-cagc behavior of glassy systems on 
a lattice and produce blocked nonergodic states and dynamical heterogeneities [18, 19]. 

Despite their success, such simple models are criticized because they have no energy relaxation, possess no underlying 
crystal phase, and fail to exhibit the correct decay of dynamical correlations with time. Experimental studies of glass 
transitions [1] are primarily presented with temperature as the control parameter, and the heat capacity plays a 
central role [2, 20]. While continuum calculations reflect many of these aspects of the system rather well [17, 21, 22], 
computer simulations of such accurate models are enormously challenging to carry out and to interpret in terms of 
simple theories. 

The aim of this paper is to bridge the gap between these two poles of scientific study by presenting a treatment 
that has the simplicity of the lattice models but schematically behaves as a true glass. The outcome of our work is a 
quite realistic model of the glass which produces a synthesis of the elements of the glass phenomenon. 

In Section II the model is introduced and the meaning of its definition is illustrated, in Section III the main 
equilibrium properties of the model are presented, in Section IV and V the dynamics and the aging behavior are 
analyzed. Finally, in Section VI the conclusions of this work are presented. 



II. THE MODEL 

One of the aims of the model is to describe crystallization with as simple a theoretical tool as possible. Continuum 
models have had extraordinary success in this area. Remarkably, simulations of hard spheres show the full range of 
phenomena which have been observed in experiments [21, 23, 24]. Here we try to represent, in a lattice model, the 
behavior of a system of quasi-hard spheres in the continuum. Our model consists of only two ingredients: a repulsive 
interaction described by a simple Hamiltonian and a kinetic rule affecting the probability of a movement. 
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A. The Repulsive Hamiltonian 



Several Hamiltonian lattice models have been introduced to mimic on a lattice the relevant properties of hard spheres 
with short ranged interactions [11, 25-29]. The equilibrium properties of our model are inspired to the Biroli-Mezard 
model [11]. The Biroli-Mezard model is defined by a many body short ranged repulsion between nearest neighbors. 
In their paper, Biroli and Mczard are mainly interested in studying the glassy behavior. Therefore, they consider a 
mixture of two types of particle in order to hinder crystallization. In out model, instead, as the dynamical slowing 
is governed by the kinetic rule, we use the Biroli-Mezard approach only to guarantee the presence of an underlying 
liquid-crystal transition at the equilibrium. Therefore, for our purposes we consider the single-type particle case. 
Moreover, in order to represent softness, we use a modified definition of the Hamiltonian, which reads: 



Here cr is the maximum number of nearest neighbors that may surround a particle without it incurring an energy 
cost, rij the number of nearest neighbors of the particle at the j-th site, Vr is the strength of the repulsive interaction, 
9(x) is the Heaviside function and V the total number of sites (volume). Of course, after this extension the model is 
no longer a-thermal. In fact, at T — the hard repulsion is recovered as a particular case. 

The meaning of this Hamiltonian is conceptually illustrated in Figure 1. If we consider an assembly of spheres in 
the continuum, the short-ranged interaction between them only starts to become important when the neighborhood 
of a considered particle is crowded. In a lattice, this concept is usually quantified by counting the number of nearest 
neighbors of a particle on a lattice site. Therefore, the threshold in the number of neighbors distinguishes between 
crowded and not-crowded environments. The softness provided by this Hamiltonian, then, in principle allows a particle 
to be completely surrounded. Thus, the model represents soft particles that can be highly packed with an energy 
cost. This cost is proportional to the repulsive potential and is higher the lower the temperature. The bidimcnsional 
sketch in Fig. 1 shows a particle in position A which has some overlap with its neighbors. The presence of such close 
neighbors implies that the particle in position A has high energy and every movement to a position without overlaps 
(for example position B in the figure) is favourable. Thus, the kinetics of the system is characterized by a tendency 
of particles to go from more to less crowded locations. 



It is known that nonergodic systems are characterized by a peculiar behavior of particle movements. Particles move 
inside a very small space for a typical characteristic time. The reason is that the surrounding particles provide a 
cage within which the movement of the considered particle is constrained. Breaking of a cage requires a long and 
therefore improbable sequence of movements. Occasionaly, the particle may escape from the cage and make longer 
movements before being trapped in another cage. These two different behaviors are identifiable in experiments [30] 
and reproducible by continuum simulations [31]. 

The Kob- Andersen model [3] correctly reproduces this intra-cage behavior of glassy systems, as well as many other 
important signatures such as blocked nonergodic states and dynamical heterogeneities [18, 19, 32]. Being a lattice gas 
at the equilibrium, this model does not present any ordered state or thermodinamic phase transition. Our aim, then, 
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B. The Kinetic Constraint 
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FIG. 1. Schematic drawing of soft disks in two dimensions. The large contour of each particle represents an effective size, so 
that soft particles can overlap. For example, the particle in position A has higher energy than it would have in position B. 
Therefore, the movement from position A to B is energetically favourable. 

is to implement a similar kinetic rule in our model and study the interplay arising between arrest and crystallization. 
Therefore, we define a kinetic constraint as follows. A particle can only move from a site i to a nearest neighboring 
site j if the following rules are satisfied at the same time: 

(a) site j is empty; 

(b) the sum of the nearest and next-nearest neighbors of the particle in i is at most equal to a fixed parameter ck; 

(c) if the particle moves to j, the sum of its new nearest and next- nearest neighbors is at most equal to Ck ■ 

This definition is based on the idea of associating the rate of motion from the central to a neighboring empty site 
with the number of particles that surround or cage it. Thus we consider that caging by a dense set of particles in the 
immediate surroundings of the central particle will last longer (because an opening will present itself less frequently) , 
and there is likely to be a threshold beneath which there are so few particles that there is no effective caging at all. 
As an example, the sketch in Fig. 2(a) represents a particle which is caged by several neighbors in such a way that the 
illustrated movement is only possible if the original cage of particles fluctuates and opens sufficiently to provide an 
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exit path. The relative unlikelihood of this originates in the fact that there are relatively few such configurations. The 
local free energy will therefore reflect this by having a barrier between these two adjacent local minima (illustrated in 
Fig. 2(b)). By definition, then, this model aims to represent the a-relaxation of glass-forming systems but does not 
describe /3-relaxation. 

The choice to also consider the next-nearest neighbors in the definition of the kinetic rule, unlike the Kob- Andersen 
model, deserves a comment. This choice is based on which particles on the lattice are caging. If we consider the 
picture of the caging mechanism in Fig. 2(a), the number of particles constituting the local cage is larger than the 
number of particles that control the amount of space available for local (intra-cage) motion of the central particle. 
The reason is that the typical size of a cage is even smaller than the radius of the particles [30], involving only the 
closest neighbors, whereas the barrier can affect a larger number of particles, since a sequential chain of movements 
might be necessary in order to break a cage. It is therefore natural that, in applying the kinetic constraint, we count 
the nearest- and next-nearest neighbors of the caged particle. This definition is also interesting for a more technical 
reason. Because of the discrete nature of the model, a wider range of values for ck allows a more satisfactory fine 
tuning. The constraint can also be made soft, giving a finite probability for unfavourable moves. 




(a) (b) 



FIG. 2. (a) Schematic representation of the caging phenomenon. To go from xa to xb the particle has to overcome the barrier 
(b) due to particles in the immediate vicinity (such as C and D). Cage escape rates, that generally depend on the type of 
interaction and particle density, are represented by kinetic rates in Kob- Andersen models [3]. A_E is the difference in local free 
energy between cages. 

This model, therefore, merges characteristics of both lattice gas models, such as the Biroli-Mczard, and kinetic 
models, as the Kob- Andersen. A recent paper presents a new Biroli-Mczard-type model with a three species mixture 
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[33]. The composition is specifically designed to avoid crystallization and exhibits properties of a fragile glass- 
forming liquid. The model successfully reproduces a stretched exponential law of the time relaxation, dynamical 
heterogeneities, Stokes-Einstein violation and other signatures. The focus of this paper is to explore a model which 
explicitly includes ingredients motivated by current physical intuitions about glass-forming liquids: namely a short- 
ranged repulsion which determines an excluded volume and a kinetic rule which is a direct model of the cage effect. 
Moreover, our model contains a genuine equilibrium crystal phase and allows a quite realistic comparison with colloidal 
systems which crystallize in the presence of low polydispcrsity [34] . Indeed, an investigation of the interplay between 
dynamical arrest and crystallization will be presented in a following publication [35] . 



C. Simulation Rules 



The two characteristics of the model can be implemented together in a Monte Carlo scheme in the following way. 
The master equation for a dynamical process on a lattice reads 

d t P(A, t) = J2 [W B ^ A P(B, t) - W A ^ B P(A, t)] , (2) 

B 

where P(A,t) is the probability that the system is in the state A at time t [36]. In our case the states A and B are 
simply sets of occupancy numbers. Wa^b is the rate of transitions from A to B. At the equilibrium, the left side of 
(2) is zero, then a condition that satisfies the master equation is that for every B: 

Pl q _ W A ^B _ p(E A -E B ) ( o\ 

P B q ~ W b ^a ~ { ) 

where E A , E B are the energies of the states A and B, respectively. The general expression for Wa^b is the product 
of a kinetic term K and an energy term F: 

W a ^b = K(A 1 B)F a ^b, (4) 

where 

( e -P(E B -E A ) ( Eb _ Ea ) >0 

Fa^b = < r (5) 

[ 1 (E B - E a ) < 

and K(A, B) is a symmetric function with respect to A and B. It is evident that the choice of the function K(A, B) 
does not affect the equilibrium state, but only the zones of the phase space visited during the simulation. We only 
consider movements involving one particle at a time, with one lattice step displacement. Thus, the only difference 
between state A and B consists of a single particle movement to an adjacent site. Without any loss of generality, we 
can label the states A and B with the site indexes i and j of the places involved in the movement. Hence, we rewrite 
the transition rate probability for a particle going from site i to site j as: 

Wi^j=K(i,j)Fi^j. (6) 

In a quite general way, we define the kinetic term as: 

K(i,j) = cxp{-V R -[(m l - c K )Q{m l - c K ) + 

+{m j - c K )6{mj - c K )]}- (7) 
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Here m, and rrij are the sums of nearest and next-nearest neighbors of the particle before and after the movement, 
respectively, 0(x) is the Heavisidc function, Vk is the barrier height (Sec Figure 2). However, the results presented in 
this work deal with the hard limit of the constraint, which can be explicitly defined as: 

K(i,j) = 6(c K -rni)6(c K -rrij), (8) 

where we assume the convention 6(0) = 0. The soft kinetic rule (7) simply provides a smoother behavior than the 
hard rule (8). Results from simulations show that the soft rule may be used to interpolate between the discrete values 
of ck , but it gives no additional interest to the model. 

The procedure described here aims to represent faithfully the kinetics of the model. The Monte Carlo scheme 
chosen is not supposed to reach the equilibrium state quickly and efficiently (in fact, it is very slow). On the contrary, 
our purpose is to represent the movements as they actually occur in a real system. The main assumption that we 
make is that a single particle scheme is able to accomplish this task, the concern being the fact that real dynamics 
involves simultaneous movements of particles. Indeed, it is known that dynamics can be realistically represented by 
physical-move kinetic Monte Carlo simulations when the degrees of freedom considered are a slow subset of all of the 
degrees of freedom [37]. This is correct in our case, especially because we consider high density systems, where the 
intra-cage movements are rare events. In lattice models, a number of lattice sweeps is able to rebuild the average 
dynamics of the system, so that the macroscopic outcome is the same. 

III. PHASE DIAGRAM AND GLASS TRANSITION 

The equilibrium phase diagram is not affected by the kinetic rules discussed above and has been already studied in 
detail [38]. Here we give a summary of its properties. For cr = 3, in a cubic lattice, the phase diagram is shown in Fig. 
3. The crystal phase is characterized by double layered diagonal planes with a periodicity of \/3 lattice steps. In the 
region involving a density roughly between 0.5 and 0.7 we observe the typical sequence of fluid, coexistence between 
fluid and crystal, and crystal. At higher density the crystal re-melts into a fluid through a first-order phase transition 
which is driven by the high number of defects in the crystalline structure, making the disordered high-density fluid 
more entropically favoured [38]. 

In this paper we will focus on cr = 3 and ck = 10. These values are chosen in order to examine the most interesting 
cases for our purposes. For other parameter values, the model contains a rich variety of phenomena that we do not 
address here. 

To illustrate effects commonly seen in real glasses, we have simulated the behavior of the temperature dependence of 
the mean energy, for different cooling rates, at constant density. For convenience's sake, we rescale the temperature to 
an adimensional variable T = kT'/Vft, and refer to this adimensional temperature in the rest of the paper. As shown 
in Fig. 4, for p = 0.69 and T > 0.9 (in the fluid phase according to Fig. 3), the system is only slightly slowed. Between 
T ~ 0.9 and T ~ 0.6, just after crossing the liquid-crystal first order transition, the energies decrease progressively and 
the system behaves as an under-cooled liquid. Below T ~ 0.6, the energy is constant and far from the crystal energy. 
Decreasing the cooling rate, the picture remains much the same, though the energy plateau becomes progressively 
lower, as expected. The states reached at T < 0.4 consist of particles that are almost completely blocked, and 
using these states as initial conditions in simulations even after long times it ~ 10 6 MCS) the configurations remain 
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FIG. 3. Section of the equilibrium phase diagram of the model for c_r = 3 [38], -F^fiuid, C=crystal. In the crystal phase for 
T < 0.4 a glassy state is observed (for ck = 10). 



In Fig. 4 we also show hysteresis between a cooling and a heating curve (at rate R = 2.510 -7 ). In a cooling process 
the relaxation time increases until we reach a certain point where it becomes longer than the cooling step. This causes 
the system to fall out of equilibrium and it freezes into an arrested state at low temperature. As the direction of the 
structural relaxation process is always towards the equilibrium, when the temperature is increased, the E — T curve 
follows a different path from the cooling curve, joining the equilibrium liquid curve after a small delay. 

The heat capacities shown in the inset of Fig. 4 are calculated from the derivatives of the energy plot. A small 
qucnch-rate dependent peak is observed. 



We study the dynamics of the model on a cubic lattice, using the van Hove self-correlation function, i.e. the 
probability that a particle has traveled distance r in time t. In a lattice, the definition of G s is: 



where 6 is the Kronccker operator. 

The spatial Fourier transform F s (k,t) of G s (r,t) is named self intermediate scattering function and contitutes a 
useful quantity to study the relaxation properties of a system. The onset of dynamical slowing is classically associated 



unchanged. 



IV. DYNAMICAL STEADY STATE BEHAVIOR 
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FIG. 4. Energy per particle with decreasing temperature. Monte Carlo simulations of a cubic lattice of size 30 with ck = 10 
and p = 0.69; lines are a guide for the eye. Time is measured in Monte Carlo Sweeps (MCS). For different cooling rates 
R = AT/ At (2.5 ■ 10~ 7 (O), 1.25 ■ 10~ 7 (□), 5 • 10~ 8 (x) ), we observe behavior typical of glass-forming liquids (see text). 
No arrest is present for ck = 18 (R = 5 • 10~ 7 , A) and the system quickly crystallizes. In one case (R = 2.5 • 10~', +) 
the temperature has been increased from the arrested state. Hysteresis is observed. The inset shows the heat capacities for 
decreasing temperature at three of the studied rates. 



with F s (k,t) being well fitted by the Kohlrausch- Williams- Watts (KWW) stretched exponential: 

/ t \ 0(T) ' 



(10) 



As one approaches the glass transition, typically r diverges and j3 decreases [1]. We calculate F s and, following a 
typical experimental approach, we fit the temperature dependence of the self intermediate scattering function by the 
KWW law. It is known that for k — > 0, F s satisfies the gaussian approximation 

1 



F s (k,t) = exp 



2d 



k 2 (Ar{ty 



(11) 



because after moving a very large distance, particle movements always become uncorrelated [2]. On the other hand, 
for k ~ > oo, i.e. at very short distances, movements are always correlated, even in a perfectly diffusing liquid. In order 
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to study the glass transition, then, we have to choose an intermediate value of k for which it is easy to check the 
validity of the KWW law (10), over the whole temperature range of interest. Of course, in a finite lattice model the 
values of k are both bounded and discrete, so that k can be neither too small nor too large. However, some values 
at small k can still be affected by the gaussian approximation. In order to check that, we consider a fluid state at 
high temperature (T = 1.5) and then we run simulations at the same density but at the temperatures studied in 
Fig. 4. In Fig. 5 plots of F s (k,t) versus k 2 are presented for fixed time t = 2 • 10 5 . In our plots the label k is the 
wavenumber scaled in units of 2-7T, so that 1/k directly corresponds to a displacement. The linear initial behavior at 
small k is compatible with the expectations. At high T, the relaxation is obviously quicker, so that the crossover to 
the high k regime is recognizable. The fact that at low T the behavior is gaussian at every k does not imply that the 
system is diffusive, because, as we will see further, the mean square displacement is not linear in time. Figs. 6 and 
7 represent the behavior of F s (k,t) versus time at fixed k = 0.577 ~ (corresponding to the crystal periodicity) 

and k ~ 0.091 (corresponding to 11 lattice spaces), respectively. From Fig. 6 it transpires that at high temperatures 
the value of F s at t = 2 • 10 5 becomes extremely small, so that the observed behavior is affected by numerical error. 
On the contrary, at k = 0.091 (Fig. 7), F s is finite at all temperatures and therefore allows an easy check of the 
stretched exponential law (10). 



r J 



O 




T=0.60 
T=0.64 
T=0.68 

T=0.70 
T=0.72 
T=0.74 
T=0.78 



T=0.85 
T=0.95 
T=1.00 

T T=1.10 
T=1.20 
T 



-+- 

- □ 

- x- 
-o- 

-A- 



-A-- 

- <$» - 

- +- 
-□- 

1.45 - -x- 



0.0 



0.2 



0.4 



0.6 



0.8 



FIG. 5. Analysis of the k dependence of F 3 (k,t) for p — 0.69; lines are a guide for the eye. An equilibrated fluid at T = 1.5 is 
used as initial state of simulations running at many temperatures from the range T = 0.60 — 1.45. All the data refer to fixed 
time t = 200000. 



The dynamical analysis we carry out is based on the assumption that the states involved are stationary. This is an 
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FIG. 6. Analysis of the time (t) dependence of F a (k,t) for fixed k — 0.5774. As in Fig. 5, temperatures are in the range 
T = 0.60 -1.45. 



approximation because the system is actually aging slowly. We check the validity of this assumption by plotting the 
energy evolution with time in the case of cooling at a fixed rate R = 2.5 ■ 10~ 6 , for a few temperatures (Fig. 8). It 
can be seen that the system is stationary for most temperatures with only a small amount of relaxation observed for 
T = 0.55. We will examine the aging properties of the model in the following section. 

In Fig. 9 the dependence of /3(T) on temperature is shown for p = 0.69. The set of data refers to a simulation 
of cooling performed taking as initial state the final state of the previous temperature. We find that significant 
non-exponential slowing arises around T — 1.0 and, consistent with expectations, the system develops non-Arrhenius 
behavior thereafter (see also inset to Fig. 10) [1, 21]. This outcome is interesting, because it illustrates the fact that 
dynamical slowing can develop long before a true glass transition. 

In order to estimate the Kauzmann temperature of the model, we assume a classical dependence of the characteristic 
time on temperature [39] : 

^roexp(^A-). (12) 

The infinite temperature relaxation time tq is given by fitting the values t(T) at 2.0 < T < 3.0 to an Arrhenius form 
[21]. In Fig. 10 the fit of the law in (12) is calculated for p = 0.69. The extrapolation of the data leads to a divergence 
at Tk = 0.42 ± 0.03 which would then be considered the glass transition. 

The extrapolated value of the Kauzmann temperature Tk corresponds in our model to a temperature where 
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FIG. 7. Analysis of the time (t) dependence of F a (k,t) for fixed k — 0.0914. As in Fig. 5, temperatures are in the range 
T = 0.60 - 1.45. 



disordered configurations remain disordered for a time scale much longer than the one of the fast processes. Only 
slow aging remains, as we are going to show further (Fig. 12), and the dynamics is sub-diffusive with no tendency 
to change (as shown by Fig. 16). This is compatible with a scenario of strong space correlations due to caging and 
vanishing configurational entropy. 

Other functional forms can fit the data with some approximation. In Fig. 11 we plot our data in the Bassler form 
t = Aexp(E/RT 2 ) [40]: the fit is quite good at high T, but it slightly worsens in approaching the supercooled region, 
as the dynamics slows down and it moves away from a normal fluid. 

Let us focus now on the mobility of the system in the glassy region of the phase diagram. We consider a fluid state 
at high temperature (T = 1.5) and then we run simulations at the same density but at much lower temperature (for 
example T — 0.4), where the crystal is the equilibrium state and at the onset of the hypothetical glass transition 
(Fig. 3). Fig. 12 shows the evolution of energies for two representative densities, on the left and on the right hand 
side of the hypothetical arrest transition. 

It is important to realize that, though only cr and Vr are relevant to the equilibrium phase behavior, all three 
parameters, including Ck, arise from underlying microscopic interactions, and a general description of maxima, saddles 
and minima of a potential surface is impossible without them. Thus, in the absence of the caging kinetic rules (ck = 18) 
the system quite quickly crystallizes (Fig. 12). Furthermore, for relatively low densities (e.g. p = 0.64) the system 
crystallizes almost immediately because caging is ineffective. On the other hand, for p = 0.69 and at low temperature 
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Energy versus time for cooling rate R = 2.5 • 10 -6 . Data from Fig. 4. 



the system evolves so slowly that crystallization has not been observed in any accessible time for our simulations 
(t ~ 10 8 ). To our knowledge, this dramatic glass-like freezing of the system is present in lattice models only when 
both short-ranged repulsion and barrier crossing effects arc combined as in the present model. 

The kinetic rule is a fundamental part of the model, as the Hamiltonian in itself does not incorporate the caging 
phenomenon and therefore the dynamical heterogeneities typical of quasi arrested systems. We believe that this 
kinetic rule is well justified by experiments which show the presence of caging and dynamical heterogeneities [30, 41], 
and our model should be considered part of this tradition. Moreover, there is extensive evidence that these types of 
kinetically constrained models reproduce a remarkable set of dynamical heterogeneities [18, 42]. 

To better capture the onset of dynamical slowing, some authors perform an inherent structure (IS) analysis which 
has been successful in both continuum and lattice models [21, 43]. However, this method appears to be suitable only 
for models which are fully described by a Hamiltonian. Following the IS implementation for lattice models illustrated 
in [43], we consider an algorithm which is going through every particle in typewriter order and try in sequence a 
movement in all the six possible directions: 

• If AE < and the kinetic rule allows it, do the move. 

• If AE = and the kinetic rule allows it, do the move with probability 1/2. 

• In any other case, reject the move. 



However, this procedure leads to a crystalline state in our model, as is shown by the IS calculation of samples at 
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FIG. 9. Stretched exponent j3 from the fitting of the KWW law (for F a (k,t)) against T, for cooling simulations at the rate 
R — 2.5 ■ 1CP 6 ; p — 0.69 and k = 0.0914. For high T the decay is simply exponential (/? ~ 1), whereas at T ~ 1 we observe a 
progressive deviation from this behavior. 



different temperatures in Fig. 13. The reason is that in our model the interplay between short-ranged repulsion 
and dynamical heterogeneities is partially destroyed by this analysis, because the systematic search of lower energy 
movements promotes unlikely paths, and destroys long-lived islands of quasi-blocked particles. In other words, the 
systematic search of the IS analysis allows us to find the very few movements which unlock the whole system but 
which would be very unlikely in a stochastic approach. 



V. AGING 



The time evolution of the energy as in Fig. 12 clearly shows that the system ages, even at very high density. 
Indeed, we know that the model presents a crystal phase at the equilibrium and therefore aging phenomena are 
expected. We have already seen that for high ck the system crystallizes. However, at ck = 10, we do not observe 
crystallization at high density. In spite of that, the system can be heterogeneous, containing little crystallites in an 
amorphous background. In order to understand better this issue, we study the static structure factor S(k). In Figures 
14 and 15 we show the evolution of the structure factor during the process of cooling. At low density (p = 0.64), the 
system quickly develops a single peak corresponding to the crystal periodicity of y/i lattice spacings. At high density 
(p = 0.69), the structure factor does not present the previous typical behavior. The profile of S(k) is much lower, as 
expected from direct observation of the configurations, but, on the other hand, it is evident that some degree of order 
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FIG. 10. Kauzmann plot for k — 0.0914. The Kauzmann temperature Tk = 0.42 ± 0.03 is obtained by extrapolation. 
Representative errors are marked on the graph, taken from tens of independent runs. Inset: Arrhenius plot (see e.g. [21]). The 
quantity plotted, Tln(r/ro), is constant when r displays Arrhenius behavior. The infinite temperature relaxation time tq is 
obtained by fitting r(T) values for T = 2.0 — 3.0 to an Arrhenius form. Deviation from the constant, high-temperature value 
is seen around T ~ 1. 

is present. 

In Fig. 16, the mean square displacement is plotted versus time for a system at p = 0.69. Up to t ~ 10 7 , the data are 
well fitted by a power law (r 2 ) ~ t 7 , with 7 = 0.344. This value is remarkably similar to the one found in experiments 
involving different types of particles [44, 45] . The fact that the agreement is also quantitatively acceptable could be 
the signature of a new universality class in the area of arrested matter. To our knowledge, this is a new achievement 
for a lattice model. The second regime presents a value 7 = 0.265 which is related to aging of the sample. Comparing 
this plot with the energy evolution in Fig. 12, it appears that the breakdown of the subdiffusive law i 34 is connected 
with a change of convexity in the energy evolution at about t w 10 7 . This regime not only shows that the system 
slows down with time, but also indicates that it evolves to a disordered, perhaps completely blocked, non-crystalline 
state. It is also possible to note that for higher densities the characteristic time separating the two power laws shifts to 
lower values. Therefore, the 7 = 0.34 law seems to be characteristic of the arrest transition, because moving beyond 
the arrest transition restricts its domain of validity. 

We can interpret the slowing of the system beyond arrest as caused by the existence of a new kind of state 
dynamically similar to random close packing. Particle movements do not lead to configurations that are more and 
more similar to the crystalline equilibrium state. On the contrary, the pathway makes the system evolve into a sort 
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FIG. 11. Bassler law fitted for low (a) and high (b) temperatures. 
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FIG. 12. Energy evolution of the model with time (ck = 10, size 30 ). Energy difference with respect to the crystal energy is 
plotted for densities p = 0.64 and p — 0.69, at T = 0.4. Starting from a fluid state, the time scale is dramatically different in 
the two cases: for p — 0.64, the system crystallizes in less than 10 5 MCS, for p = 0.69, an extremely slow evolution is observed. 
Setting ck = 18 for p — 0.69, i.e. removing the kinetic constraint, the system very quickly crystallizes. 



of side-track: the slowing process is progressive, perhaps leading to a completely arrested state. 

Our interpretation of a subdiffusive exponent close to 1/3, based on the analysis above, is that the phenomenon 
happens only when the free energy landscape has some particular characteristics. First, it has to present a global 
minimum corresponding to a stable crystal phase. As a consequence, the temperature has to be low enough to 
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FIG. 13. Inherent structure energy relaxation of samples at different temperatures T. 



guarantee quasi-hard particle interaction. Second, there must be many local minima, separated by high barriers, so 
that the time evolution to the global minimum requires long times to overcome or circumvent the barriers. Therefore, 
it is clear that, in order to observe the phenomenon, the parameters of the model have to be tuned in such a way that 
these conditions are satisfied. It is still questionable how the value of the subdiffusive exponent (which could be a 
rational number) can be theoretically justified. In this paper, we only show that this law, observed in a few different 
experiments, can even be reproduced by a simple lattice model, pointing out that the phenomenon seems to be a 
robust property of the arrest transition. 

As a last remark, from Figure 16 it emerges that the mean square displacement can be relatively high. After 10 8 
time steps, we have yj (r 2 ) ss 2.2, which means that on average every particle has traveled more than 2 lattice steps. 
In spite of setting T Tk, this should not be surprising, because a glassy behavior is perfectly compatible with 
subdiffusivity. The fact that the diffusion constant is zero does not mean that every particle is blocked, but only that 
there is not a characteristic time in the distribution of wait times per particle [46] . 



VI. CONCLUSIONS 



The model may be used to study many other interesting properties, but discussion of the details are not central 
to our point here. Rather we note that, for the first time, a simple lattice model has been able to reproduce the 
range of phenomena from real glasses and energy landscape models, including rate dependent slowing of the energy, 
divergence of the characteristic time with an appropriate KWW law, vanishing of the diffusion constant, onset of 
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FIG. 14. Static structure factor versus time for a cooling simulation at p = 0.64 and rate R — 2.5 • 10 6 
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FIG. 15. Static structure factor versus time for the cooling simulation at p = 0.69 and rate R = 2.5 ■ 10 as in Fig. 4. 
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collective behavior, and a possible quantitative representation of subdiffusivity in dense systems. It is interesting to 
explore the origins of these new effects in the model. 

Models based on kinetic constraints alone create a complex effective free energy landscape in which, at sufficiently 
high density, many movements are prohibited by infinite or large barriers. Nonetheless, many dynamical pathways 
involving long ranged transport still remain at fixed (zero) energy. True dynamical arrest only occurs in such models 
when the lattice is fully filled. This dramatic reduction of dynamical pathways induced by kinetic constraints certainly 
leads to dynamical slowing, but not to true glassy behavior. If we now allow different local energies within different 
cages one obtains a complex energy landscape. Then, rare pathways that were formerly barrier-less remain 'easy', 
but acquire a multitude of smaller energy barriers. The accumulation of such bumps against the backdrop of a 
vanishing number of easy pathways ultimately leads to interesting singular behavior for the characteristic time, which 
is considered to be truly representative of the glassy state. This is the root of glass behavior, and its physical origins 
are quite clear in our model. 

Regarding the behavior of the mean square displacement, lattice models clearly show that sub-diffusion arises when 
particle movements are intermittent and there is not a characteristic wait time [46]. In a dense system, subsequent 
movements are allowed by the same hole, coming back according to a random walk scheme. Therefore, arbitrary long 
return paths cause the absence of a characteristic time. This behavior changes when particles are being moved by 
several different holes, and we observe diffusion, or crystallization occurs. However, this behavior does not break down 
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at the arrest transition. If the system has diverging subdiffusive domains, there is not a time scale for their breaking. 
Moreover, the kinetic properties of our model are due to the underlining Kob- Andersen kinetic rules, which have been 
thoroughly studied in the past [3]. In such kind of kinetic models, it is possible to identify a dynamical correlation 
length which describes the typical size of domains where local movements of particles do not contribute to diffusion 
[19]. When this length diverges, the system becomes dynamically arrested. The presence of a characteristic length 
in our model could give an explanation for the remarkable behavior of colloids under gravity, where the gravitational 
length does not capture the nature of the phenomenon [44] . 

It must be considered intriguing that the ingredients for such a model are already known in the literature but 
that they have not hitherto been combined in this way. Purely repulsive Hamiltonian lattice models without kinetical 
barriers [11, 25, 47] appear not to yield a KWW characteristic time law, as in experiments and continuous simulations. 
On the other hand, purely kinetic models do not have a crystal phase [3]. Here we have a simple model that reproduces 
the main effects associated with glassy systems. Thereby it opens up the possibility of a more transparent dialogue 
between experimental scientists, simulators in off-lattice models, and theorists working in such highly simplified 
models. 

We acknowledge with pleasure discussions with W. Kegel and S. Granick. This paper is written within an EU-US 
research consortium, funded in part by a grant from the Marie Curie program of the European Union, ('Arrested 
Matter', contract number MRTN-CT-2003-504712), and by an SFI grant. 
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